Effect of particle inertia on the turbulence in a suspension 
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We propose a one-fluid analytical model for a turbulently flowing dilute suspension, based on 
a modified Navier-Stokes equation with a fc-dependent effective density of suspension, p c s(k), and 
an additional damping term oc 7 P (fc), representing the fluid-particle friction (described by Stokes 
law). The statistical description of turbulence within the model is simplified by a modification 
of the usual closure procedure based on the Richardson-Kolmogorov picture of turbulence with a 
differential approximation for the energy transfer term. The resulting ordinary differential equation 
for the energy budget is solved analytically for various important limiting cases and numerically in 
the general case. In the inertial interval of scales we describe analytically two competing effects: the 
' energy suppression due to the fluid particle friction and the energy enhancement during the cascade 

process due to decrease of the effective density of the small scale motions. An additional suppression 
or enhancement of the energy density may occur in the viscous subrange, caused by the variation 
of the extent of the inertial interval due to the combined effect of the fluid-particle friction and the 
decrease of the kinematic viscosity of the suspensions. The analytical description of the complicated 
interplay of these effects supported by numerical calculations is presented. Our findings allow one 
to rationalize the qualitative picture of the isotropic homogeneous turbulence of dilute suspensions 
as observed in direct numerical simulations. 
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Introduction 



The interaction of solid particles or liquid droplets with the turbulence in a gas controls the 
performance of various engineering devices and is important for many practical applications like 
the combustion of pulverized coal and liquid sprays and cyclone separation. This interaction plays 
also an important role in many areas of environmental science and physics of the atmosphere. Dust 
storms, rain triggering, dusting and spraying for agricultural or forestry purposes, preparation and 
processing of aerosols are typical examples. For a review of turbulent flows with particles and 
droplets see, e.g. the book by C.T. Crowe, M. Sommerfeld and Y.Tsuji [1]. 

In dilute suspensions with small volume fractions of particles, C p , the particle-particle interactions 
are negligible. Nevertheless, for p p j pf ^> 1 (the ratio of the solid particle material and the gas den- 
sities), the mass loading <f> = C p p p / p{ may exceed unity and the kinetic energies of the particles and 
the carrier gas may be comarable. Hence the "two-way coupling" effect of the fluid on the particles 
and vice versa must be accounted for. Current understanding of the turbulence in dilute suspensions 
is still at its infancy due to the highly nonlinear nature of the physically relevant interactions and a 
wide spectrum of governing parameters (the particle size a vs. L and rj, the outer and inner scales of 
turbulence, the particle response time r p vs. j L and j v , the turnover frequencies of L- and rj- scale 
eddies). 

Existing analytical studies of the problem are mainly based upon a two-fluid model description 
wherein both the carrying fluid and particle phases are treated as interpenetrating continua [1-4]. 
This model deals with non-interacting solid spherical particles with a radius a small enough such 
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that: 

1. One can neglect the effect of preferential concentration and may assume homogeneity of the 
particle space distribution. This is not always so. Above some critical radius a cr the space 
homogeneous distribution of particles becomes unstable. Resulting clustering instability leads 
to preferential concentration. For a detailed theory of this effect, see Ref. [5] and references 
therein. In the present paper we consider only particles with a < a cr . 

2. The Stokes viscous drag law for particle acceleration, du p /dt = [ut — w p ]/r p , is valid (itf is the 
fluid velocity). 

Unfortunately, the statistical description of two-fluid turbulence with closure procedures requires a 
set of additional questionable simplifications due to the lack of understanding of the relevant physics 
of the particle-fluid interactions. This makes closures of the two-fluid model highly qualitative at 
best [4, 6, 7]. 

We think that the basic physics of the problem may be better described by a simpler one-fluid 
model for turbulent dilute suspensions, which uses standard closure relations of one-phase turbulence. 
The present paper suggests such a model and, as a first step, uses a properly modified simple 
closure, based on the Kolmogorov- Richardson cascade picture of turbulence. The resulting non-linear 
differential equation for the energy budget were solved analytically. This provides an economical 
and internally consistent analytical description of the turbulence modification by particles including 
the dependence of suppression or enhancement of the turbulence on the three governing parameters: 
(r P 7 L ), <f) and the scale of eddies. These effects were previously observed in numerous experimental 
and numerical publications, see, e.g. the review by Crowe, Trout and Chung [8]. Many groups 
carried out experimental work; for an overview see Pietryga [9]. Other researchers studied the 
modification of turbulence by small particles using direct numerical simulations [4, 10-13] or by 
large-eddy simulation [14]. Nevertheless the complicated physics of turbulently flowing suspensions 
in the two-way coupling regime still wait for a detailed analytical description. 

Our analytical findings in this paper successfully correlate important features of turbulence mod- 
ification observed in numerical simulations Refs. [4, 12, 13]. We believe that the one-fluid model 
(together with more advanced closures of one-phase turbulence) offers an insight in basic physics of 
particle-laden turbulent flows. The next step in this development should include the effect of pref- 
erential concentrations, which was studied so far only for a given turbulent flow field of the carrier 
fluid [5]. 

The paper is organized as follows. In Sec. II we review, after a presentation of the notation 
(Sec. I A) and an evaluation of the characteristic time scales (Sec. IB), some publications about 
DNS-simulations (Sec. II A), about experimental work (Sec. II B) and about some analytical mod- 
els(Sec. II C). A critical evaluation of the existing analytical models[4, 15-20] is made. 

In Sec. Ill we suggest a new one-fluid equation of motion (3.1) for turbulently flowing suspensions 
with small particles. This is a modified version of the Navier-Stokes equation with two wave-number- 
dependent parameters, p e s(k) and 7 P (fc): 

• The /c-dependent effective density of suspensions p e s(k) describes the different degree of involve- 
ment of heavy particles in turbulent fluctuations with different wave-numbers [referred to below 
as fc-eddies]. For /c-eddies with a turnover time l/^(k), which is much smaller than the particle 
response time r p , the particles may be considered at rest and p e s(k) is about the density of the 
fluid itself, p f . For /c-eddies with T p ^y(k) <C 1 the effect of the particle inertia may be neglected 
and particles may be considered as fully involved in the motion of eddies. Therefore for small 
enough k the effective density p e s{k) is close to the mean density of the suspension (fluid plus 
particles), p s = pf(l + 0). Our Eq. (3.2) reasonably describes p e s(k) for all values of k. 
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• The damping term 7 P (/c), given by Eq. (3.3) describes the fluid-particle viscous friction. The 
function j p (k) saturates at the level l/r p for small scale eddies with r^ik) ^> 1, when the 
particles may be considered to be almost at rest. In this regime the damping is /c-independent, 
while the turnover frequency of /c-eddies j(k) grows with k. Therefore for large k j p (k) -C j(k) 
and the particle-induced damping of these /c-eddies may be neglected with respect to their 
energy loss in the cascade process, which is determined by the frequency 7(fc). In contrast, for 
small enough k [when r p ^{k) <C 1] the particles are almost completely involved in the motions of 
/c-eddies and their contribution to 7 P (/c) is suppressed by the factor [r P 7(A:)] 2 <C 1 with respect 
to 1/V P . 

Our one-fluid model for turbulent suspensions (3.1) is first postulated in Sec. Ill A. Its physical 
interpretation is discussed in Sec. Ill B. A detailed derivation of Eq. (3.1) is given in Sees. Ill C, 
HID and HIE. The most difficult problem here is how to account for the nonlinear effect of the 
interaction of /c-eddies within the one-fluid model of turbulent suspensions. The suggested form of 
the nonlinear term (3.5) is a modification of the standard Navier-Stokes nonlinearity and is based 
on: 

• a rigorous description of eddy interactions in both limiting cases T p ^{k) <C 1 and T p ^{k) ^> 1 

• respect of the fundamental symmetries of the problem - Galilean invariance and conservation 
of energy. 

Section IV deals with the budget of the kinetic energy in turbulently flowing suspensions. One has 
to account not only for the dissipation of energy due to the fluid-particle friction but also for the 
effect of particles on the energy redistribution in the system due to the eddy interaction. First we 
derive in IV A the budget equation (4.1) which accounts for the energy pumping due to a stirring 
force, energy damping due to the kinematic viscosity and fluid-particle friction and also describes 
the flux of energy over the scales due to the nonlinearity of the problem. Equation (4.1) is exact 
but unfortunately is not closed. As usual it includes a 3rd order velocity correlation functions. As a 
first step in the analysis of turbulent suspensions in the framework of our one-fluid model Eq. (3.1) 
and the budget Eq. (4.1), we use in this paper, sect. IV B, a simple closure procedure based on 
the Richardson-Kolmogorov cascade picture of turbulence in which the energy flux is accounted for 
in a differential approximation. Needless to say that there are various closure procedures for the 
Navier-Stokes turbulence in the literature. They may be straightforwardly applied to our Eq. (3.1). 
This important part of the project will be done elsewhere. 

The derived energy balance equations are summarized in Sec. IV C. They have a very simple and 
transparent analytical form (4.22) - (4.26), allowing their effective analytical analysis, see Sects. V 
and VI. In particular in Sec. VB we found a simple solution for the case of micro-particles having 
a very small response time. In Sec. V C we found the iterative solution for the case of a suspension 
with heavy particles in the inertial interval and analyzed its accuracy in Sec. VD. 

In section VI we analytically describe a complicated interplay between two competitive effects: of 
the turbulence suppression and the turbulence enhancement in the inertial interval of scales, as well 
as in the viscous subrange. A brief comparison of our finding with DNS results is done in Sec. VI E. 

In the concluding Sec. VII we summarize the results of the paper and present our ideas for further 
work. 
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I. NOTATIONS AND RELEVANT TIMESCALES 
A. Nomenclature 

• p{, u(t, r) - density and velocity of the fluid 

• u(u>,r), u(t,k), u{uj,k) - Fourier transform of u(t,r) with respect to t [t — > u], with respect 
to r [r — > k], and to both variables [t — > u , r — > k] 

• F(t, fe), F(u, k) - pair correlation functions of fluid velocity in (t, k) and (u, k) representation 

• E(k) = p{k 2 F(0, k)/2n - one dimensional spectrum of the turbulent kinetic energy of the pure 
fluid (fluid without particles) 

• £(k) - one dimensional spectrum (of the turbulent kinetic energy) of the suspension 

• j(k) - turnover frequency of fc-eddies (turbulent fluctuations of the characteristic scale l/k). 
May be understood also as l/r(k), where r(k) is the life time of /c-eddies. In the Kolmogorov 41 
picture of turbulence 'y(k) ~ k^kE(k)/p { . 

• E = J dkE(k)/2n, £ = J dk£(k)/2n - total turbulent kinetic energy of respectively the pure 
fluid and the suspension 

• a, p p , m p = 47ra 3 p p /3 - radius, density and mass of the particles 

• C p , £ 3 = 1/Cp , volume fraction of particles and volume of suspension per particle 

• ■?/>= [4na 3 /3}/£ 3 , <fi = m p /p { £ 3 - volume fraction and mass loading parameter 

• r p - particle response time, also referred to as Stokes time scale 

• t l - turnover time of the energy containing eddies (of scale L) 

• 5 = t p /t l - the particle response time in the units of r L . 

= rj/u. n - Kolmogorov (viscous) microscale; characteristic velocity and time at scale rj 
of turbulence 

• Pcs(k) ~ effective density of the suspension for turbulent fluctuations of characteristic scale l/k 
[referred to as /c-eddies] 

• v, v e e(k) - kinematic viscosity of the pure fluid, effective kinematic viscosity of /c-eddies in the 
suspension 

• 7 P (A;) - effective damping frequency in the suspension due to the fluid-particle friction 

• e(k) - (one dimensional) flux of the turbulent kinetic energy of the suspension via a sphere of 
radius k in fc-space, also referred to as energy flux over scales. 

B. Evaluation of time scales 

The radius of the particles is supposed to be small enough, so that the particle Reynolds number 
7\fc p is less than a critical value (73e cr ). In this case we can apply the Stokes approximation (according 
to which the fluid-particle friction force is proportional to the difference between the particle velocity 
and the fluid velocity). Careful analysis by Lumley [21] shows that in a turbulent flow the condition 
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for the validity of 7& p ^ 7Ze CT may be expressed via the particle radius a and the Kolmogorov 
micro-scale r\ in the following way 

a < 2r?(p f /p p ) 1 / 3 . (1.1) 

It is clear that one of the important parameters in the physics of turbulently flowing suspensions 
is the ratio of the inertial time scale of the particles (the Stokes time scale r p ) and the life time t v 
of eddies of the Kolmogorov micro-scale. The particle response time is given by 

r P = a = ' y 1 - 2 ) 

07T v p{ a 9 pf v 

where we use the expression for the particle mass m p : 

4:71 

As is well-known the Kolmogorov micro-scale f] is found from the condition that the Reynolds number 
for eddies of scale r\ is equal to unity: 

-Re 11 = V v v /u=l . (1.4) 

Here v v is the characteristic velocity of 77-scale eddies. It is related to the turnover time of these 
eddies in the following manner t v = r)/v v . This allows us to rewrite the requirement (1.4) as follows 

Tr, = V 2 /"- (1-5) 

The ratio of the time-scales r p and t v immediately follows from Eqs. (1.2) and (1.5): 



Zk - - ?E — 

r v 9 pi rf ' 

Substituting the condition (1.1) for the validity of the Stokes approximation we find 

/ \ i/3 
Z£ < I ^£ 



(1.6) 



(1.7) 



where we neglected the difference between 8/9 and 1. Equation (1.7) means that for "heavy" particles 
in a gas, that satisfy Stokes approximation, the particle response time scale may be about ten times 
larger than the Kolmogorov time scale: r p < 10r v . For such particles in a liquid the two time scales 
are about the same. So we may conclude that for "heavy" particles in a gas, that satisfy Stokes 
approximation, the inertia of the particles may be expected to be important in a considerable part 
of the energy spectrum. For particles in a liquid the particle inertia will only be significant for the 
smallest eddies, for which r p ~ t v . 



II. REVIEW OF LITERATURE 



This section is devoted to a review of the literature about the problem of a turbulently flowing 
suspension. We will review important findings from published numerical experiments, physical 
experiments and analytical models. 
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A. Review of some DNS simulations 

To study the two-way coupling effect several groups have applied the direct-numerical-simulation 
technique (DNS) to particle-laden isotropic turbulence. A brief review of some of the publications 
is given below. 

Squires and Eaton [10] considered the particle motion in the Stokes regime in which gravitational 
settling was neglected. They assumed statistically stationary isotropic turbulence. Mass loadings 
from zero to unity were considered for a series of particle response times varying from 0.3 t T] to 11 t v , 
where t v is the Kolmogorov time scale. They found that the overall reduction in turbulence kinetic 
energy for increasing mass loading was insensitive to the particle response time. They attributed the 
non-uniform distortion of the turbulent energy spectrum by particles to the preferential concentration 
of particles into regions of low vorticity and/or high strain rate. 

Elghobashi and Truesdell [11] examined turbulence modulation by particles in decaying isotropic 
turbulence. They used the particle equation of motion derived by Maxey and Riley [23], and found 
that for the large density ratio considered in their simulations the particle motion was influenced 
mostly by drag and gravity. They found that the coupling between particles and fluid resulted in 
an increase in small-scale energy. This increase in the energy of the high-wave-number components 
of the velocity field resulted in a larger dissipation rate. They also found that the effect of gravity 
resulted in an anisotropic modulation of the turbulence and an enhancement of turbulence energy 
levels in the direction aligned with gravity. 

Boivin, Simonin and Squires [4] also made a very detailed DNS-study of the modulation of isotropic 
turbulence by particles. The focus of their work was on the class of dilute flows in which particle 
volume fractions and inter-particle collisions are negligible. Gravitational settling was also neglected 
and the particle motion was assumed to be governed by drag with particle response times ranging 
from the Kolmogorov scale to the Eulerian time scale of the turbulence and particle mass loadings up 
to unity. The velocity field was made statistically stationary by forcing the low wave-numbers of the 
flow. Like in [10, 11] the effect of particles on the turbulence was included by using the point-force 
approximation. The DNS-results showed that particles increasingly dissipate fluid kinetic energy 
with increased mass loading, with the reduction in kinetic energy being relatively independent of the 
particle response time (as was already found in [10]). The viscous dissipation in the fluid decreases 
with increased mass loading and is larger for particles with smaller response times. The fluid energy 
spectra show that there is a non-uniform distortion of the turbulence spectrum with a relative 
increase in small-scale energy (as was found in [11]). They state that the fluid drags the particles at 
low wave-numbers while the converse is true at high wave-numbers for small particles. 

Sundaram and Collins [12] performed DNS-simulations of particle-laden isotropic decaying turbu- 
lence. The particle response time was in the range 1.6 r rj < r, p < 6.4 t v . The ratio of the particle 
density and fluid density was of the order 10 3 . The particle Reynolds number 73e p remained less than 
7?e cr , and the drag force on the particles was described by Stokes law. The point-force approximation 
was employed to represent the two-way coupling force in the fluid momentum equation. The DNS- 
results showed that the particles reduce the turbulent kinetic energy as compared to the particle-free 
case, and this reduction is less pronounced for smaller response times t p . The results also showed 
that the total turbulent energy dissipation is increased by the particles, and the increase is larger for 
smaller t p . The turbulent energy spectrum is reduced at small wave-numbers and increased at high 
wave-numbers by the two-way coupling, and the location of the cross-over point is shifted towards 
larger wave-numbers for larger r v . 

Druzhinin [13] examined the modulation of isotropic decaying turbulence by microparticles, for 
which 2a < rj, r p < t v and 7^e p < 73e cr . The gravitational settling is neglected. Due to the fact that 
p p ^> pf, the mass loading may be large enough to modify the carrier flow. Druzhinin first derived an 
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approximate analytical solution for the energy spectrum and then performed also DNS simulations. 
The results obtained for particles whose r p < 0At v show that both the turbulence kinetic energy 
and the turbulence dissipation rate are increased by the two-way coupling effect as compared to the 
particle-free case. For particles with sufficiently high inertia (r p > 0.5t v ) the two-way coupling effect 
caused a reduction in the turbulence kinetic energy as compared to the particle free case. Druzhinin, 
therefore, showed that there occurs a qualitative transition in the two-way coupling effect of particles 
on isotropic turbulence as the particle response time is increased from r p <C t v (microparticles) to 
r p ~ r v (particles with finite inertia). For microparticles there is an increase of all wave- numbers in 
the energy spectrum. For particles with a higher inertia that is no longer the case. 

B. Review of some laboratory experiments 

Many experiments have been carried out to study the modulation of turbulence in the carrier 
phase by particles. An overview of the experimental work up to 1999 is given by Pietryga [9]. 
Experimental measurements in shear flows, e.g. particle-laden jets and boundary layers, have shown 
that turbulence velocity fluctuations may be either increased or decreased due to the modulation 
of the flow by (heavy) particles. However in turbulent shear flows it is often difficult to separate 
the direct modulation of the turbulence due to the momentum exchange with the particles from the 
indirect changes occurring through modification of turbulence production mechanisms via interaction 
with mean gradients. In grid-generated turbulence these production mechanisms are absent. It 
approximates in the best possible way the homogeneous, isotropic turbulence with particles that we 
study in this publication. We will, therefore, briefly review below some literature publications about 
experimental work devoted to the study of the modulation by particles of grid-generated turbulence. 

Schreck and Kleis [24] studied the effect of almost neutrally buoyant plastic particles (density 1050 
kg/m 3 ) and heavy glass particles (density 2400 kg/m 3 ) on grid-generated turbulence in a water flow 
facility. The average particle size was 655yum. The particle Reynolds number of the plastic particles 
Ve p ~ 8, for the glass particles 7\fc p ~ 20. The particle volume fraction was varied between 0.4% 
and 1.5%, so the system was very dilute. Mean velocity and velocity fluctuations of both phases 
were measured by a laser-Doppler velocimeter. The presence of the particles in sufficiently high 
concentration modified the turbulence downstream of the grid. The decay rate of the turbulence 
energy increased monotonically with particle concentration. The additional dissipation rate for the 
suspensions with the heavier glass particles was about double that of the almost neutrally buoyant 
plastic particles. A simple model based on the slip velocity between the phases under-predicted the 
measured increase in the dissipation rate. Schreck and Kleis, therefore, assumed that a large portion 
of the additional dissipation is associated with the measured modification of the spectral distribution 
of the turbulence energy. They speculate that the particles enhance the transfer of energy to smaller 
eddies extending the dissipation spectrum to smaller scale. Since only part of the high wave-number 
end of the spectrum could be resolved experimentally, this speculation could not be conclusively 
demonstrated by their experimental data. 

Hussainov et al. [25] studied the modulation of grid-generated turbulence by coarse glass particles 
in a vertical downward channel flow of air. Two different types of grids were used. Glass beads with 
an average diameter of 700/im and a mass loading of 10% were used. The particles were about 7 
times larger than the Kolmogorov length scale r\ and 7& p ~ 70 or 93, dependent on the type of grid 
used. The particle response time scale of the particles r p was about 5000 to 7000 times larger than 
the Kolmogorov time scale r v . The mean velociy and the turbulence intensity along the channel axis 
(and in some cross-sections) were measured by means of a laser-Doppler velocimeter. The decay 
curves of the turbulence intensity in the streamwise direction showed an attenuation of turbulence 
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intensity of the flow by the particles. The particles caused an increase in the total dissipation rate 
of the turbulence. Hussainov et al. found that the presence of the particles decreased the energy 
spectra at high frequencies. This seems to be in contradiction with the speculation of Schreck and 
Kleis, that the particles enhance the transfer of energy to smaller eddies. 



C. Analytical models 

The starting point for analytical models, described in the literature, is often the Navier-Stokes 
(NS) equation for the velocity of the pure fluid (fluid without particles) u(t, r) 



Pi 



d +(u- V)-uV 2 }u + Vp = f p + f, (2.1) 



idt 



where p(t, r) is the pressure and pf is the fluid density. The random vector field f(t, r) represents 
the stirring force responsible for the maintenance of the turbulent flow. Equation (2.1) includes also 
the force f p (t, r) caused by the friction of the fluid with particles: 

f p (t,r) = ^[v(t,r)-u(t,r)} . (2.2) 

Here v(t,r) is the velocity field of the particles, considered as a continuous medium with density 
m p /£ 3 = pf<fi, where m p is the mass of a particle, £ 3 - suspension volume per particle and <f> is the 
mass loading parameter 

<P = m p / Pi e. (2.3) 

The validity to represent f p (t, r) in the form of (2.2) is based on the assumption of space homogeneity 
of the particle distribution. It is also assumed that the particles are small enough for the Stokes 
drag law to be valid. The equation of motion, suggested in the literature, for the continuum phase 
of the particles does often not include the pressure and viscous terms 



| + (v ■ V) v = -f p . (2.4) 



Equations (2.1) and (2.4) were used by Baw and Peskin [15] to derive a set of "energy balance" 
equations for the following functions: 

• Eff(k) - energy spectrum of the fluid turbulence, E(k) in our nomenclature 

• Efi iP (k) - energy spectrum of the fluid turbulence along a particle trajectory 

• E{ p (k) - fluid-particle covariance spectrum 

• E pp {k) - particle energy spectrum 

In the balance equations the following energy transfer functions occur 

• T ff f(/c) - energy transfer in fluid turbulence 

• T iP: f(k) - transfer of fluid-particle correlated motion by the fluid turbulence along the particle 
path 

• T iP:P (k) - transfer of fluid-particle correlated motion by the particles 

• T pp {k) - transfer of particle-particle correlated motion by the particle motion 
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• II qi f(A;) - fluid-particle energy exchange rate. 

Baw and Peskin [15] made a set of simplifying assumptions in order to be able to analyze the balance 
equations. First, they assumed that the particles do not respond to the fluid velocity fluctuations 
due to their (very large) inertia. Therefore 

Eg, p (k) = E s (k) , (2.5) 

Ttp, f (k) = T lP:P (k) = Tpp jP (k) = . 

This assumption is, of course, not realistic for particles satisfying the Stokes' approximation. Their 
next assumption 

n q j = <!>[E fp (k)-E B , v (k)]/T p , (2.6) 

may be understood as a statement that the fluid-particle exchange rate is statistically the same for 
all scales characterized by a /c-independent frequency 7 P = 0/r p . This is reasonable for particles with 
very large inertia, but then Stokes law is not valid. For particles satisfying Stokes law, assumption 
(2.6) has to be replaced with a more realistic, /c-dependent frequency 7 P (&)- We will come back to 
this point when discussing our new model. 

A serious difficulty in the derivation of the energy balance equations is how to find a closure 
expression for third-order velocity correlation functions, responsible for the various energy transfer 
functions. Baw and Peskin assumed that Tg${k) can be expressed similarly as in the case of a pure 
(single phase) flow 

d e) /3 k 5 / 3 E s (k) 

where 6f is the viscous dissipation in the pure fluid (without particles) and a is the so-called 
Kolmogorov constant. This assumption seems questionable to us. According to the spirit of the 
Richardson-Kolmogorov cascade picture of turbulence one may express inertial range objects, like 
Tg t f(k) in terms of again inertial range quantities, like k, Eg{k) (which is done in Eq. (2.7)) and 
e(k), the energy flux in fc-space. In a single phase flow, indeed e(k) = 6f. However this is not the 
case for a turbulent suspension due to the fluid-particle energy exchange, given by Eq. (2.6). We 
think that our closure (to be discussed later on) is an improvement in this respect. 

With this simplified model Baw and Peskin predicted the following influences on the energy spec- 
trum of the fluid turbulence due to the particles: 

• a decrease of the energy in the energy-containing range of the spectrum 

• an increase in the inertial range of the spectrum 

• a decrease in the viscous dissipation range. 

Boivin, Simonin and Squires [4] used the same model as in Ref. 15. They also applied assumptions 
similar to Eq. (2.6) and Eq. (2.7). Fortunately, they took into account the response of the particles 
to the turbulent velocity fluctuations by relaxing assumptions (2.5) and also accounted for the very 
important physical effect of the energy dissipation due to the drag around the particles. For that 
reason they approximated Tg t f(k) and Tf P) f(fc) as follows: 

T s>i (k) = --^ [et - U q>i (k)} 1/3 k 5 / 3 E s (k) , (BSS1) 
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T{ ^ {k) = -^djb^-^*)] 178 * 573 ^^) • 

Notice that this closure has the same weakness as Eq. (2.7), involving the dissipation range value ef. 
With the above described changes with respect to the model described in Ref. 15 Boivin, Simonin 
and Squires found an increase in the viscous dissipation range of the fluid turbulence spectrum for 
small values of the particle response time r p . 

Al Taweel [16] calculated the rate of additional energy dissipation due to the presence of the 
particles. Because of their inertia the particles were assumed not to follow completely the turbulent 
velocity fluctuations of the carrier fluid. They expressed the additional dissipation in terms of the 
turbulent kinetic energy of the suspension. Then they added this term to the balance equation of 
the turbulent kinetic energy, making the (questionable) assumption that the energy flux across the 
spectrum has the same functional form as in a single-phase flow. Solving this equation they found 
an attenuation of the high-frequency fluctuations with a small alteration of the energy-containing 
low frequencies. Although there was an additional energy dissipation due to the particles, the total 
energy dissipation was reduced due to the reduction of viscous dissipation in the carrier fluid. 

In a number of publications [17-20] Felderhof, Ooms and Jansen developed an analytical model 
for the dynamics of a suspension of solid spherical particles in an incompressible fluid based on the 
linearized version of the Navier-Stokes equation. In particular they studied the effect of the particles- 
fluid interaction on the effective transport coefficients and on the turbulent energy spectrum of the 
suspension. Also the hydrodynamic interaction between the particles and the influence of the finite 
size of the particles were incorporated. However it is needless to say that the nonlinearity of the 
Navier-Stokes equation is of crucial importance in the problem of turbulence. Felderhof, Ooms and 
Jansen were well aware of this problem, but wanted to study in particular the influences of the 
particle-particle hydrodynamic interaction and of the finite particle size at a high particle volume 
concentration. 



III. ONE-FLUID MODEL NAVIER-STOKES EQUATION FOR TURBULENT SUSPENSIONS 

In Sec. II C we discussed the two-fluid model of suspensions consisting of the Navier-Stokes (NS) 
equation (2.1) for the fluid and Eq. (2.4) for the "gaseous" phase of particles. This approximation 
is based on the assumptions of space homogeneity of the particle distribution and applicability of 
the Stokes drag law for the fluid-particle friction. We think that the basic physics of a turbulently 
flowing suspension with these assumptions may be described in the framework of the much more 
simple one-fluid equation. This model is presented in Sec. Ill A, discussed in Sec. Ill B and "derived" 
in Sees. HID and HIE. 



A. The model 



The following equation may be considered as a model equation for turbulently flowing suspensions: 

Pcs(k) + 7 p(A0 + 7 o(fc)]«(*, fe) (3.1) 
= -Af{u, u} t ,k + f(t, k) , 
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The linear part of this equation involves: 



,„„;(//) = fH {l-t + 4> * + 2 ™W }. (3.2) 

[1 + r p7 (A;)J 

7pl J (l + ^)[l + 2r p7 (fc)] + [r P 7(fc)] 2 ' 1 } 

lo (k) = u cS (k)k 2 , u eS {k) = . (3.4) 



PcffO) 

The nonlinear term in Eq. (3.1) has the usual NS equation form: 



(3.5) 



However the vertex r^ fc differs from the standard vertex r y^ lk2 of the NS equation (see 
e.g. Refs. 26, 27): 

722k = ^[P a "(k)k^ + P^(k)k^5(k + k 1 + k 2 ), 

(3.6) 

as follows: 

r a/3 7 _ „ f 2 hk 2 k 3 \ 7fc£fc 2 { o 7 \ 

Our model differs from the standard NS equation in three aspects: 

a. Equation Eq. (3.1) involves the /c-dependent effective density of suspensions p c s{k) given by 
Eq. (3.2). The function p c s(k) satisfies the inequality pf < p c s(k) < pf(l + 0). One could say 
that p e s(k) — pf represents the contribution of the particles involved in turbulent fluctuations with 
characteristic scale 1/k to the effective density of suspensions. 

b. Equation (3.1) includes the additional damping term 7 P (fc), Eq. (3.3), describing the loss of 
kinetic energy caused by the viscous fluid-particle friction. 

c. In the absence of a stirring force f(t, r) and both damping terms, Eq. (3.1) conserves the total 
kinetic energy of suspensions £ [given by Eq. (3.42)] which is different from the kinetic energy E of 
the fluid itself. 

The explicit form (3.5) of the nonlinear term is not necessary for the simple closure procedure that 
we applied in this publication. For the introduction of the energy flux in used closure procedure it 
is enough to use the fact that the modelled nonlinearity must be conservative. However, the explicit 
form is needed for more advanced closure procedures that we intend to use in future work. For this 
reason we include it in this publication. 



B. Physical interpretation of the one-fluid model 

In a simplified fashion we may interpret p e s(k), the /c-dependent density of suspension in our model 
equation (3.1) as follows. 

Denote as / com (fc) the fraction of particles co-moving with the fc-eddies (turbulent fluctuation with 
some wave-number k) , in the sense that their velocity is the same as the velocity of /c-eddies. These 
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particles also participate in the motion of eddies with smaller wave-number k' < k but not necessarily 
in the motion of /c"-eddies with k" > k. For small k the turnover frequency 7 (fc) of /c-eddies is small 
in the sense 7(fc)r p <C 1. Therefore, in this region of k, the particle velocity is very close to that of 
the carrier fluid and we can describe the suspension as a single fluid with effective density p c s(k), 
which is close to the density of suspension: 

Ps = p f (l - V) + C pPp = p f (l - i> + 0) , (3.8) 
iP = C p [Ana 3 /3], = C p p p /p f . 

Here C p is the particle concentration, ip and are the volume fraction and mass loading parameter. 
However, for large k, when 7(A:)t p 1, the particles cannot follow these very fast motions and may 
be considered at rest. Thus the particles do not contribute to the effective density and p e s(k) — > pf. 
In the general case p e s(k) may be written as 

PcS {k) =p f [l-^ + 0/com(fc)], (3.9) 

Here a statistical ensemble of all particles, partially involved in the motion of A;-eddies, is replaced 
by two sub-ensembles of "fully co-moving" (fraction f com (k) ) and "fully at rest" (fraction f TCSt (k) = 
1 _ /com(fc)) particles, which does not contribute to p e s(k). 

The particles at rest cause the fluid-particle friction. According to Newton's third law, the damping 
frequency of a suspension 7 P (fc) may be related to the particle response time, r p , via the ratio of 
total mass of particles M p at rest to the total effective mass of the suspension M cS (k): 

m M P C p p p f rest (k) _ (f) p { f rcs t(k) 

7pl ' r p M cS (k) r pPcS (k) r pPcS (k) ■ 1 • ) 

As we mentioned, the fractions f com (k) and f res t{k) depend on r P 7(fc). Moreover, the portion f rest (k) 
is independent on the sign of the velocity, therefore we expect f Tes t(k) ~ ( r p7(^)] 2 - 111 the opposite 
case, when l/r P 7(A;) is small, / com (k) has corresponding smallness: f CO m(k) ~ l/r P 7(fc). As a simple 
model of such a function we adopt 

AestW = 1 " /com(A:) = [r p7 (A;)/(l + r p7 (A;))] 2 . (3.11) 

Using Eq. (3.11), we rewrite Eqs. (3.9) and (3.10) as Eqs. (3.2) and (3.3). Note that these equations, 
which follow from the physical reasoning described above, give the same expression for 7 P (A;) as 
Eq. (3.3) in our "derivation" in Sec. HIE 2. We consider this fact as a strong support of the 
physical relevance of our one-fluid model for a turbulently flowing suspension given by Eqs. (3.1)- 
(3.4), with /c-dependent effective density, fluid-particle damping frequency 7p , and effective kinematic 
viscosity v e &{k). 



C. Basic assumptions 

The theory developed in this paper is based on a number of assumptions and simplifications 
described below: 

1. All particles in the suspension are spheres with the same density p p and the same radius a. 

2. The radius of the particles is small enough, see Eq(l.l). 

3. The particle-particle interaction will be neglected, assuming that the volume fraction »/; C 1. 
Nevertheless, for the very heavy particles with p p ^> pf, the mass loading may be of the order 
of unity, leading to a significant modification of the turbulent flow by particles. 
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4. The turbulent flow is stationary, homogeneous and isotropic. 

5. In our equations for the energy balance (4.1) we will use simple (but physically relevant) closure 
procedures based on our effective (one-fluid) Navier-Stokes equation for suspensions (3.1) and 
on the Richardson-Kolmogorov cascade picture of turbulence. 



D. Formal derivation of the effective NS equation for suspensions 



In the derivation we begin with the NS Eq. (2.1) for the fluid component. Instead of the averaged 
expression (2.2) for the fluid-particle friction force we will use the following detailed expression 



Mt,r)=Y,* , P (t*r j )6(r-r j ), 



(3.12) 



in which F p (rj,t) is the force between the fluid and j-particle positioned at r — Tj. Assume [as in 
derivation of Eqs. (2.2) and (2.4)] that the statistics of particles is independent of the statistics of 
turbulence and, moreover, that their distribution is space homogeneous. In that case, we can replace 
the sum over the position of particles by a space integration: 



E 



1 



drj, 



where £ 3 is the volume per particle. In this approximation 

f p (r,t) = F p (r,t)/£ 3 . 



(3.13) 



We compute F p (t, r) for small enough particles with a radius a satisfying inequality (1.1), such that 
the fluid flow in the vicinity of a particle may be considered as laminar (assumption III C-2). Then 
one can apply Stokes law for the force F p (t, r): 



F p (t,r) = {[v p (t)-u(t,r)] 
with the friction coefficient ( for heavy particles (with the density p p ^> pf) is given by 

( = 6n pf v a . 
The Newton equation for a particle reads: 

dv p (t) 



A formal solution of this equation 



= -F p (t,r) = C[u(t,r)-v p (t)} . 



v P (t) 



u(t,r) , 



allows one to express the force F p (t, r) as follows 



F P (t,r) 

Here r p is the particle response time: 



m 



T 'Tt + 1 



u(t,r) . 



t p = m p /6ii v pi a 



(3.14) 
(3.15) 

(3.16) 
(3.17) 

(3.18) 
(3.19) 
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The total time derivative (d / dt) as usual takes into account the time dependence of the particle 
coordinate r: 



d 
dt 



(3.20) 



Due to the (immersed) particle inertia they do not follow Lagrangian trajectories of fluid particles. 
Therefore, generally speaking, (d/dt) does not coincide with the total Lagrangian time derivative 
in the fluid, 



D 
~Dt 



d_ 

dt 



+ u(t,r) ■ V 



(3.21) 



Consider the relationship between (d /dt) and (D/Dt): 

Du(t,r) du(t,r 



Dt 
du(t, r) 

dt 
d 1 



dt 



+ iv«-u a (t,r)]V a u(t,r) 



-u(t,rrV a lUl 
dn + r PdI 



dt 1 



Pdt 



d_ 

'dt 



l + r p —)+C 



u(t,r) , 



t u(t, r) = t p [(v Pt i ■ V)u(t, r) - (u ■ Vi)«i] . 

Here u\ = u(ti,ri), Vi = d/dri, and all derivatives with respect ti and 7*1 are taken at t\ 
ri = r. Together with Eq. (3.21) this gives: 

Du(t,r) 1 



Dt 



l + r p | + £ 



u(t,r) 



(3.22) 

(3.23) 
: t and 

(3.24) 



For particles with a small response time Ferry and Balachandar [22] show, that the particle velocity 
depends only on local fluid quantities (the velocity and its spatial and temporal derivatives). They 
derive an expansion of the particle velocity in terms of the particle response time which generalizes 
those of previous researchers. For large values of the ratio of the particle density and the fluid density 
and for small values of the particle response time our Eq. (3.24) for the force on a particle gives the 
same equation for the particle velocity as derived in Ref . [22] . 
Substitution of Eq. (3.24) into NS equation (2.1) yields: 



Pf 



d 

dt + 



- V) ]{ 1+ rr^)}^^ 

= Pi uV 2 u + f, 



(3.25) 



where is the mass loading parameter. For simplicity we consider here only the case of heavy parti- 
cles with negligibly small volume loading parameter, ^ < 1. However, the mass loading parameter 
may be of the order of unity. For example, for the water droplets in the air, p p / pf m 10 3 and for 
0=1, the volume fraction if) pa 1CT 3 . 

The inverse operator in Eq. (3.25) may be understood as a Taylor expansion with respect to the 
nonlinearity (u ■ V): 



C 



l + r p (| + £) l + r p | (1 + Tp-l)- 



+ . . . 



(3.26) 
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This expansion produces higher order [ in (it- V) ] nonlinear terms in the effective NS equation (3.25). 
These terms are not important for big eddies with r P 7(A;) <C 1 for which the operator in the braces 
in the LHS of Eq. (3.25), {...}, is close to the factor 1 + 0. In the opposite case, for small scale 
eddies with T p ^(k) 3> 1 the operator {. . . } = 1. Both limiting cases one easily gets from the first 
term in the Taylor expansion (3.26) in which there is no contribution from C. It means that only for 
intermediate scales with r P 7(A;) ~ 1 this operator may be quantitatively important. For a qualitative 
description of the "transient" process between these two regimes it is enough to account for the first 
term of expansion (3.26). In this approximation the turbulent fluid velocity around the particle is 
approximated by the velocity at a fixed coordinate, which is reasonable in statistical sense and exact 
in the limit r p ^(k) <C 1. With this approximation Eq. (3.25) turns into 



Pcflf 



d i 

+ («• V)Ju + Vp = Pi vV 2 u + f , (3.27) 
1 + -^}, (3.28) 



dt 

Pcflf = Pf 



where p c s may be considered as an operator of effective density for a suspensions. 

Since we are interested in the incompressible flows, we can project the potential components out 
of the equation of motion. This may be done by the projection operator "P, defined via its kernel 

V^{r) = J -^P a(S (k) exp[-ifc ■ r] , (3.29) 
P a(3 (k) = 5 a/3 - k a k^/k 2 . (3.30) 
The application of "P to any given vector field a(r) is nonlocal, and it has the form 

[V ■ a(r)] a = J dr{P a \r - n)a^(n) . (3.31) 



Applying V to Eq. (3.27) we find 

d 



PcflF 



r)t +V-(u-V) u = P( uV 2 u + f, (3.32) 



This equation together with the definition (3.28) for the operator of the effective density constitutes a 
one-fluid description of a turbulently flowing suspension. However, the operator form of the effective 
density is not convenient for practical calculations. To overcome this inconvenience we will derive 
below another form for the effective parameters of this equation. 



E. NS equation for suspensions with fc-dependent parameters 



In our analytical description of space homogeneous, stationary turbulence it is convenient to con- 
sider Eq. (3.32) in the (fc, u;)-representation: 

MPeffH] -ipf[^ 2 + %(u)]}u(u,k) (3.33) 
= -Af{u, w} w , fc + f(w, k) . 
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Here 



/ 



u(u, k) 
Pcfr(^) = Re{p cfr (cj)} = pt 



di dr u(i, r) exp(i cj t + ik ■ r) , 



1 + 







Pf 



= Pf 



Im [{PcffM]} = j 





l + (c^r p )2j' 

(j)U 2 T p 



+ Mn 



(3.34) 
(3.35) 
(3.36) 



1 + 



1 — i UJT n 



f\f{u,u}u :k = [p cS (uj)V ■ (u • V)«] fc 



(3.37) 



A/*{w, ii} w ,fc denote the nonlinear term in (cu, ^-representation and frequency uk 2 describes the 
viscous damping. 

The Navier-Stokes equation for suspensions (3.33) involves a frequency-dependent effective density 
of suspensions, p' cS {uj) and ^-dependent frequency %(ou) responsible for the damping due to fluid- 
particle friction. To use standard closure procedures in the statistical description of turbulence 
one needs frequency independent coefficients in the basic equation of motion. On other hand, these 
procedures may be applied to equations with /c-dependent coefficients. Therefore for further analysis 
it is much more convenient to deal with a A;-dependent effective density p c s{k) of /c-eddies. To relate 
these functions we note that the /c-eddies have a characteristic frequency of motions, j(k) [related 
to their life-time r(k) by a simple relation ^y(k) ~ l/r(k)]. 



1. k-dependent effective density of suspensions 

In the inertial interval of scales Eq. (3.33) must preserve the total kinetic energy of a suspension 
£ if one neglects the fluid-particle friction 7 P (^) — > 0. The equation for £ may be written in terms 
of the density [p' eS {uj) and F(u, k), the pair correlation function of the (cu, fe)-Fourier components of 
velocity, u(u,k). Namely 

*=/If£^.*). (3 ' 38) 

(2-K) 4 5(k — ki)5(u — uji)F(uj, k) = {u{u, k) ■ u{u\, fei)) . 

For isotropic turbulence F(u, k) = F(u, k) and Eq. (3.38) allows one to introduce the one- 
dimensional energy spectrum of suspension, £{k) according to: 

^S(k), (3.39) 

Define a /c-dependent effective density of suspension, which gives the same one- dimensional spec- 
trum £{k) as the a;-dependent effective density p' eS (u): 

pMk) = SUfM^ . (3 . 41) 

J F(u, k) au 
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Then, Eq. (3.40) takes the form 



Notice that 



/ 



■^F{u,k)=F{k), (3.43) 

is the simultaneous velocity pair correlation function. With this notations Eq. (3.42) may be written 

as 

£(k) = P^k 2 F{k) , (3.44) 
while the traditional notation for 1- dimensional spectrum of kinetic energy of fluid itself is E(k): 

E(k) = ^k 2 F(k) . (3.45) 

2,7V 

Formally speaking, in order to evaluate p e s(k) by Eq. (3.41) we need to know the lo dependence 
of F(u,k). This is not a simple task. Instead we will use a few reasonable forms of F(ui,k) and 
compare the resulting functions p e s(k). One of the frequently used is the Lorentzian form 

F(u,k) = F(k) ^{iw ( 3 - 46 ) 

which corresponds to the simplest "one-pole" approximation for the Greens functions. Using this 
tu-dependence in Eq. (3.41) we have the following simple form for p e s(k). 

<"* W = 4 + Tt4*)]' (3 ' 47) 

For small r p j(k) this gives a correction linear in T p j(k) 

p cS {k) «p f [l-0r P 7(fc)] , (3.48) 

which contradicts the physical intuition. Indeed, one may consider /c-eddies as having a motion 
with the characteristic frequency j(k) and expect that p e s(k) may be obtained from p' eS {u) with the 
substitution to — > •y(k). This gives a correction quadratic in T p ^{k) 



p cS (k) - Pf (l + 0) » -0p f r p V(A;) • (3.49) 

This contradiction follows from the fact that the model function Eq. (3.46) decays very slowly for 
lo — > oo, like 1/lo 2 . It is known from the diagrammatic analysis of the different time velocity 
correlation function F(r,k) that for small r the difference F(r,k) — F(0,k) does not contains |r| 
and decays not slower than r 2 . Therefore the Fourier transform of F(r,k), F(u,k) has to decay 
with LO faster than 1/lo 2 , at least as 1/cj 4 . To meet this requirement we consider instead of Eq. (3.46) 
the function: 

F(u, k) = F(k) 273(fc)/7r , (3.50) 

> 2 + 7 2 (A0] 2 
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which gives instead of Eq. (3.47) 

Peff(^) = PfU + — — 72- J (3.51) 

"> [1 + r P 7(A;)J > 



p f (i + 4>) - (f>pt 



vr( fc ) v 

1 + r p7 (/c)J 



Now the correction to p e s{k) is quadratic in r p r y(k) which agrees with the expectation (3.49). One 
observes the same agreement for any other model dependence F(u, k) decaying even faster than 
1/u 4 . 

Therefore the linear part of Eq. (3.33) may be modelled as 

{up eS (k) - \pi\vk 2 + 7 P (w)] }u(uj, k) = ... , (3.52) 
with p e s(k) given by Eq. (3.51). 



2. Effective fluid-particle damping frequency 7 P (fc) 

Using Eq. (3.33) or Eq. (3.52) together with Eq. (3.40) we can compute the contribution of the 
fluid-particle friction to the damping of £{k): 

d£(k) f du 7pH, 2 



P -^)W)^r kFM - (3 - 53) 



dt 

Introduce an ^-independent fluid-particle damping friction by a standard relation: 

d£(k) 



= -2 lv {k)£{k) . (3.54) 

p 



dt 

Combining these two equations with Eq. (3.40) one gets 

= Pi f%(u)F(u,k)du 

J P'es{u)F{u,k)du 

Substitution of Eqs. (3.36) and (3.46) into Eq. (3.55) gives Eq. (3.3) for 7 P (fc). With this knowledge, 
(3.52) may be further simplified as follows 

p eS (k){uj - i [v eS (k) k 2 + 7 P (fc)] }u(lu, k) (3.56) 
= -M{u, u} Wife + / . 

Here v e s{k) is given by Eq. (3.4). Notice that this equation gives the same dissipation rate (3.54) 
due to fluid-particle friction as Eq. (3.33) and the same dissipation rate 



d£{k) 



dt 



-2u cS k 2 £(k) (3.57) 



due to the kinematic viscosity. 
The suggested form of M{u, u}^^ in terms of p e s(k) will be discussed in the following Sec. Ill E 3. 
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3. w -independent nonlinearity of the effective NS equation 

a. Nonlinearity in the usual NS equation. Consider first the nonlinear term in the "usual" NS equation 
for single-phase flow. In (u, ^-representation it has the form (see e.g. Refs. 26, 27): 

^{u, u}Z :k = J (3.58) 

X 7fciifc 2 M K^l' fc l) S 7( W 2,fc2)5(^ + CJ 1 + U) 2 ) . 

Here ^ k ^ klk2 is the so called vertex of interaction given by Eq. (3.6). It includes transversal projectors 
accounting for the incompressibility of the fluid, delta function of k- vectors originating from the space 
homogeneity of the problem and is proportional to k (as a reflection of V operator in the nonlinear 
term in r-representation). 

The vertex 7^ fc2 satisfies so-called Jacobi identity 

7fcfclfc 2 + ^fcafcfc! + lklk 2 k = > (3.59) 

as a consequence of the energy conservation by the Euler equation. 

b. Nonlinearity in the effective NS Eq. (3.1). A rigorous derivation of the nonlinear term in the effective 
NS Eq. (3.1) is a very delicate issue. For example, in Eq. (3.27) we used the operator of the effective 
density (3.28) containing only the first term of expansion (3.26). This approximation does not 
account for all terms, second orderin u. This derivation is beyond the scope of this paper. Instead 
we present here physical arguments which allows us to propose a form of M{u, u}u,k which satisfies 
all needed requirements. 

By analogy with Eq. (3.58) we can write: 

\rt ia f d 3 k l d 3 k 2 duj 1 duj 2 
{«, u} W)fc = j — (3.60) 

x r fcfc1fc 2 M K^i ) k!)u*{uj 2 , k 2 )5{u + u 2 ) , 

where the vertex T^ fc2 differs from 7^^, Eq. (3.6), because now pf ^ p c s{k). 

The simplest possible generalization of the vertex, just a replacement pf — > p c g(k) in Eq. (3.6) 
leads to a vertex ^ kk J lk2 , which does not satisfy the Jacobi identity 

1 fcfcifc 2 ' 1 fc 2 fcfci 1 fcifc 2 fc — u ) l°- ui J 

leading to violation of the conservation of the kinetic energy of suspension 8. We suggest Eq. (3.7) 
for r^ fea . Clearly, due to Eq. (3.59) this vertex satisfies requirement (3.61) and consequently, 
Eq. (3.1) conserves the energy £. 

Another physical requirement is Galilean invariance of the problem. This is the case for the 
standard NS equation with vertex (3.6) in which pf is /c-independent. For the /c-dependent density 
in the vertex (3.7) Galilean invariance implies that in the limit, when one of the wave-numbers is 
much smaller then two others (say k\ <^ k 2 = k 3 ), the effective density must depend on the smallest 
k- vector. Obviously, this is the case for the A;-argument of p c g(k) in Eq. (3.7). This guarantees 
Galilean invariance of Eq. (3.1). 

Now Eq. (3.56) involves only ^-independent coefficients and may be rewritten in (t, k) represen- 
tation, see Eq. (3.1). 
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IV. BUDGET OF KINETIC ENERGY IN TURBULENT SUSPENSIONS 

In this section we consider the budget of kinetic energy in turbulent suspension. In Sec. IV A we 
will use the one-fluid model (3.1) to derive (for homogeneous, isotropic turbulence) the following 
budget equation for the (1-dimensional) density of kinetic energy: 

^%p-+bo(k)+7 P (k)]£(t,k) (4.1) 
= W(t, k) + J(t, k) . 

The left hand side (LHS) of this equation includes two damping terms, j (k)£(t, k), caused by the 
effective kinematic viscosity and 7 p (fc)£ (t, k) caused by the fluid-particle friction. The density S(t, k) 
is given by Eq. (4.7). The right hand side (RHS) of Eq. (4.1) includes the source of energy W(t, k), 
localized in the energy containing interval, and the energy redistribution term, J7"(t, k). 

The budget equation (4.1) is exact, but unfortunately not closed. Equations Eq. (3.4) for the 
effective kinematic viscosity and Eq. (3.3) for 7 P (fc) includes "turnover frequency" of /c-eddies, 7(/c). 
Also W(t, k) contains yet unknown (u, f) correlations, Eq. (4.6). And finally J{t,k) is given by 
Eqs. (4.9), (4.3) via 3-rd order velocity correlations F 3 . There are many reasonable closure procedures 
for the approximation of high-order velocity correlations by lower-order ones. To elucidate the basic 
physics of the problem at hand, in this paper we will use the simplest possible closure. Application 
of more sophisticated closures will be done elsewhere. 

A. The energy budget equation 

In order to derive Eq. (4.1) we multiply Eq. (3.1) by u(t, k') and average: 

p^ k ){^r+bo(k)+i P (k)] F (^ k )} 

= J(t,k) + W(t,k), (4.2) 
J (t, k) = j ^0^K^ k Ff\ t] k, kl , k 2 ) . (4.3) 

Here F(t,k), and F 3 (t; . . .) are the 2nd and 3rd order simultaneous velocity correlation functions 
taken at overall time t: 

(27i) 3 5(k + fei)F(i; k) = (u(t, k) ■ u{t, k x )) , (4.4) 
(27r) 3 5(fc + fex + k 2 )Ff\t; k, fex, k 2 ) (4.5) 
= (« a (t,fc)^(i,fe 1 )«T(i,fc 2 )) . 

Note that the time t in the argument of F(k) in Eqs. (3.43) - (3.46) is omitted: F(t, k) = F(k). 
In Eq. (4.2) we introduce also the simultaneous (u, /)-cross correlation functions 

(27r) 3 5(fc - kjWit, fe) = (u(t, k)f(t, fcO) . (4.6) 

We can rewrite Eq. (3.42) for the density of the kinetic energy of suspension in terms of F(t, k) = 
F(t,k) (for isotropic turbulence): 

£(t,k) = ^k*F(t,k). (4.7) 
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Multiplying Eq. (4.3) by k 2 /2n one gets finally the balance Eq. (4.1) in which 

W(t,k) = ^-W(t,k), (4.8) 

J(t,k) = ^J(t,k) . (4.9) 

Notice that effective vertex r^ fe in Eq. (4.3) was constructed such that the total kinetic energy 

£ is the integral of motion (neglecting pumping and damping): J °° J{t, k)dk = 0. Therefore the 
energy redistribution term J(t, k) may be written in the divergent form: 

J(t,k) = _^£M>, wll ere (4.10) 
dk 



e(t,k) 



POO 

/ dkJ(t,k) (4.11) 
Jk 



is the (one-dimensional) energy flux over scales. 

In the rest of the paper we will consider only stationary solutions of Eq. (4.1). Omitting (here and 
below) time argument one finally has: 

[ 7o (*) +7p(*)]£(*) = W(k) - ^ • (4.12) 



B. Simple closure for the energy budget equation 

The effective density in our approach, Eq. (3.2), depends on the characteristic frequency fc-eddies 
7(/c). This object may be evaluated as the inverse life time of these eddies which is determined by 
their viscous damping and energy loss in the cascade processes. Accordingly, 7 (/c) is a sum of two 
terms 

7(A;)= 70 (A;)+7 C (A;), (4.13) 

where jo(k), Eq. (3.4), is the viscous frequency and 7c (A;) may be evaluated as the turnover frequency 
of A;-eddies: 

7 c (fc) ~ kU k , where U k ~ ^kS(k)/ PcS (k) (4.14) 

is the characteristic velocity of /c-eddies. We therefore define 

lc (k) = C^k3£(k)/ Pcfi (k), (4.15) 

where C 7 is some dimensional constant, presumably of the order of unity. Clearly, the same evalua- 
tion (4.15) one gets from a dimensional reconstruction of ^(k) in terms of the only relevant (in the 
K41 picture of turbulence) variables k, £{k) and p e s(k). 

In the same manner, by the dimensional reasoning, one gets the following evaluation of the energy 
flux: 

e(k) = C e y/k5£Z(k)/p eS (k), (4.16) 

where C £ = 0(1). Notice that in pure fluids [with p e s(k) — > pf] Eqs. (4.15) and (4.16) are nothing 
but the K41 evaluation of the corresponding objects. This become even more transparent if one 
rewrites Eq. (4.16) in the more familiar form: 

S(k) = d [e\k) p cS {k)] V V 5 / 3 , Ci = Cfl\ (4.17) 
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Together with Eq. (4.15) this gives a useful evaluation of 7 c (fc) via e(k): 

7c(fc) = C 2 



£(/c) fc 2 l 1/3 C 7 
, ^2 = 



Peff(^) -I 



a 1/3 



(4.18) 



Lastly, we have to evaluate the energy input in the system W(k). It follows from Eq. (3.1) that 
u(k) may be evaluated as f(k)/^(k)p e s(k). Together with Eqs. (4.6), (4.8) and (4.15) this gives: 



W(k) = C w fi 



2 jkS(k) 



(4.19) 



where C w = 0(1) and f% is the pair correlation of the forcing [which may be defined similarly 
to (4.6)]. 



C. Dimensionless budget equation 

For the convenience of the reader we present here the full set of equations which will be studied 
below. To non-dimensionalize this equation we define a dimensionless wave-number, k, and the 
integral-scale related parameters 

1 



k = k L , e L = e(-) , j L = 7Q) , p L = Peii(j^ 



Define also the dimensionless functions 



e(fc) 

! 

Pes(k) 



In = 



7(fc) 

1 

W(k) 



(4.20) 



(4.21) 



p L W(l/L) ' 

in which the argument k is written as a subscript to distinguish them from the corresponding 
dimensional functions of the dimensional argument k. 
The resulting dimensionless budget equation reads 



^ + -CT K 
dy y 

+ £(§f (1+rJ - w " 

4> 5 7 K 

K = (1 + 0)(1 + 25 7k ) + (57 k ) 2 ' 
C = C 1 C 2 , 5 = 7 L r p . 



(4.22) 



(4.23) 



Here we used Eq. (4.17) and defined the Reynolds numbers for the carrier fluid, 7\fcf, and the effective 
Reynolds number for the suspension TZe s 



Re s = 
Pl = Pi 



Lv, 



v 

Lv, 



7) 



1/3 



(4.24) 



1 + 



Pi 

1 + 25' 



(4.25) 
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in terms of the rms turbulent velocity, v L , dominated by L-eddies. Obviously, 73ef involves the 
kinematic viscosity of the carrier fluid, u, while 7Ze s depends on the effective kinematic viscosity of 
the suspension, z/ cff (L _1 ), for the outer scale of turbulence, L. 

Eq. (4.22) has to be considered together with equations for p K and 7 K , which follows from Eqs. (3.2), 
(4.13) and (4.18): 



Pk 



Ik 




C 2 TZe s p K 



(4.26) 



These two equations allow us to express the functions 7 K and p K in terms of e K . With these solutions, 
Eq. (4.22) becomes an ordinary differential equation for the only function e K . 

The first line of Eq. (4.22) describes the effect of particles in the inertial integral of scales. This 
part involves the mass loading, 0, the dimensionless particle response time (normalized by the eddy 
life time), 5, and the parameter C, characterizing our version of the K41 closure. 

The second line of Eq. (4.22) represents the effect of the viscous friction, which is proportional to 
l/7?e s , and the pumping term W K , which we choose as follows: 



1 



2wa 



exp 



(y 



2 a 2 



(4.27) 



This function has a maximum at y — 1 (the input of energy is largest at k = l/L), while the 
parameter a describes the characteristic width of the pumping region. In addition, the function W K 
satisfies the normalization constrain 



oo 

/ 



W«dy = l, 



(4.28) 



which follows from Eq. (4.22) in the limit <j < 1. 



V. SOLUTION OF THE BUDGET EQUATION 
A. Simplification of the energy pumping term 

First notice, that the turbulence statistics in the energy containing range, k L = y ~ 1 is not 
universal and depends on the type of energy pumping, in our case, on the function W K . Therefore 
for general analysis, which is not aimed at the study of some particular type of turbulence generation, 
we can take the pumping of energy in a narrow shell in the fc-space. This means 

Urn {WJ = <*(«) , (5.1) 

a— > 

where 5(k) is the Dirac 5 function. In this limit and with zero boundary conditions for e K , 7 K at 
k = (and, consequently, p K — 1 at k — 0) , Eq. (4.22) can be solved on the interval < k < 1. 
This gives 



e K = 1 , 7 K = 1 , p K = 1 , at k = 1 



(5.2) 
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In the limit (5.1), Eq. (4.22) has zero RHS for k > 1: 



+ -C T K + ^ (^f) V3 (1 + T K ) = , (5.3) 



Relations (5.2) can be considered as the boundary conditions for Eq. (5.3) at k — 1. 



B. Particle free case and limit of small particles 

Consider now the particle-free case, = 0, and the case of very small particles, 5 = 0, for finite 
7Ze s . In both cases Eq. (5.3) becomes: 

dy Ve s K J 

We took here in account that according to Eq. (4.26) p K — 1 for = and also for 5 = 0. Notice, that 
for = 0, v L — v and, consequently, 7Ze s = T&i, while for 5 = 0, v L = vj (1 + 0) and 7Ze s = 7fef(l + 0). 
The reason is that for 5 — > all particles are fully involved in turbulent motions and one can consider 
a suspension as a single but heavier fluid with the density pf(l + 0). 
The solution of Eq. (5.4) with the boundary condition e 1 — 1 is 

In the particle free case = and this solution turns into 

where 7?e s ^> Tfef as we discussed above. In the bulk of the inertial interval these solutions give a 
small viscous correction to the K41 solution with the constant energy flux e K — 1. Namely 

e K w 1 + 3 d (1 - K 4/3 )/(4-Re s ) , for K « l/7?e s . (5.7) 

The local in the /c-space closure procedure, used in the paper, works reasonably well in the inertial 
interval, where the energy exchange between eddies is dominated by the eddies of compatible scales. 
However, it is violated in the bulk of the viscous subrange, where the dynamics of eddies of very small 
scales is dominated not by their self-interaction, but by their energy exchange with 77-eddies of the 
Kolmogorov micro-scale 77. Therefore we cannot expect Eq. (5.3) to reproduce the exponential decay 
of the energy flux in the viscous subrange. Nevertheless, this equation describes on a qualitative 
level the behavior of the energy flux until the very end of the inertial interval giving the crossover 
scale to the viscous subrange, i.e. the value of 77. According to Eq. (5.5), the energy flux become 
zero at 

k cv = (1+47^/d) 3 / 4 . (5.8) 

It convenient to introduce here an effective Reynolds number of the carrier fluid and suspensions 

-Ref = 47& f /Ci , 7&f = 47& S /Ci , (5.9) 

which enters in the corresponding Kolmogorov-41 evaluations of the viscous cutoff. For example for 
the fluid 

m = L/K cr ^L/[Reff 4 . (5.10) 
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C. Iterative solution in the inertial interval 



In the bulk of the inertial interval, after neglecting the viscous terms (i.e. for 7Ze s 
becomes 



de K 



+ 



dn K (1 + 0)(1 + 257 k ) + (5 7k )2 



0, 



oo), Eq. (5.3) 

(5.11) 
(5.12) 



PK = 



1 + 



'(i + ^) 2 



1 + 



1 + 25 



1. Large scale solution of the budget equation 

In region of large scales, k ~ 1, the functions p K ~ 1 and we can simplify Eq. (5.12) by the 
replacement p K =>- 1 in the equation for j K , i.e. 7 K =>- el/ 3 K 2 ^ 3 . In the denominator of Eq. (5.11), where 
the ^-dependence of 7 K is less essential, we can make further simplification, replacing 7 K 
The resulting equation allows separation of variables: 



K 



2/3 



2 de K 



1/3 



5 d/t 2 / 3 
*o(«) 



C*o(«), 



(5.13) 



;i + 0)(i + 2<5/t 2 / 3 ) + 5 2 /t 4 / 3 ' 



The solution of this equation with the boundary conditions £i 

1 



1 is e K = e 0jK , where 



J («) 



[l+CJ (/€)] d 



(5.14) 



In 



Sk 2 ^ + 1 + - y/gj + 0) 
5+1+0-^0(1+0) 



-In 



^ 2 / 3 + l + 0+ y/0(l + 0) 
5+1 + 0+^0(1+0) 

Now with Eq. (5.12) we find the following approximations 



Po,k 



"fO,K 



1 + 



l + 254 / K V/ 3 



(i+* e ;i s «=»/3) : 

(^0,k K VP0,k) 1/3 • 



0. 



1 + 2(5 

(i+^r 



(5.15) 



Formally speaking, the analytical solution (5.14-5.15) is valid only for To find the solution 

of the initial Eq. (5.11) in the whole interval of scales, we will iterate this equation, taking the 
analytical solution as a starting function. Comparing in Sec. VD this solution with the next order 
iterations and with the numerical solution of Eq. (5.11), we will find an actual region of applicability 
of the analytical solution (5.14-5.15). 



26 



2. First improvement and subsequent iterations 



With the analytical solution (5.14-5.15), we can improve approximation (5.13) of Eq. (5.11) by 
Pk =>• Po,k [instead of p K =>- 1], which gives 7 K =>• e}/ 3 W 2 / 3 / p , e in the numerator of Eq. (5.11). In the 
denominator we replace 7 K =>• e l J 3 K 2 ^ 3 /po,e- The improved simplification of Eq. (5.11) reads: 



2 d< 



'1/3 



5 d/t 2 / 3 
*i(«) 



= C*i(«), 







P^ 3 [(l + 0)(l + 257o, K ) + (57o,«) Z ] 



Integration of this equation gives the first iterative solution of Eq. (5.11), e K = £i ;K , where 

1 



[l + CJ^/c)] 



3 ' 



J!(k) = ^ ^!(x)da; 2/3 . 



1 + 



This allows further improvement of approximations (5.15) 

l + 2(5 7 o, K 

Pl ' K = + 77T^ — V* 

7i,« = (ei,«« 2 //»i,«) 1/3 • 
Now, the next iteration steps are obvious. The n-order solution is 

1 



1 + 25 



E-n,n 

J n (n) = 

*n(«) = 



[l + CJ n ( K )] 3 ' 

^ | *„(x)dX 2 / 3 , 







+ 0)(1 + 2*7n-M) + (*7»-l,i.) ] 

j. 1 + 26 
+ 5 

(i + «) 2 



2i ' 



1 + 



1 + 257 



n— 1,k 



(l + <J7„_i,«)' 

7ra,K i^n,K^ I Pu,k) 



(5.16) 



(5.17) 



(5.18) 



(5.19) 



(5.20) 



D. Accuracy of the iterative solutions 

To get an understanding of the accuracy of the analytical solution e 0)K , Eq. (5.14), and the first 
iterative solution e\ tK , we compare them with the "numerically exact" solutions of Eq. (5.11), e K in 
the wide inertial range of four decades. 

We found that for all values of k and 5 the analytical function Eq >k works unexpectedly well for 
C 4> < 0.25. To illustrate this, we plot in Fig. 1 functions eq jK , £i )K and e K for C = 0.25 (panel a) 
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and C — 1 (panel b) for = 0.25, 0.5 and = 1 with 5 = 0.1. The relative difference between e , K 
and 5 K is about a few percents for all three cases in panel a and for the case = 0.25 in panel b. 

For larger values of the product C the accuracy of few percents is achieved in the smaller region 
of n, where r p ^f(k) = 5k 2 ^ 3 < 1, i.e. approximately for k < o~~ 3 / 2 . For example for 5 = 0.01 this is 
three decades, k < 10 3 , while for S — 0.1 only for k < 30, as we show in Fig. lb. Moreover, the first 
iterative solution, ei jK gives a very good approximation to e K for all reasonable values of parameters. 
This is illustrated in Fig. 16 for C — 1 and = 1. Notice, that for C = 0.25 and = 1 (Fig. la) the 
plots of £i )K e K are undistinguished within the line width. 

The conclusion is that for the qualitative and semi-quantitative description of the turbulence 
modification by particles in the inertial interval we can use the analytical solution (5.14-5.15), 
corrected, if needed, by the first iteration, ei jK . 



VI. TURBULENCE MODIFICATION BY PARTICLES 



A. Preliminaries 



Consider now separately the density of kinetic energy of the carrier fluid, E{(k), and that of the 
particle, E p (k) (i.e. the density of the kinetic energy of the particle velocity field). According to 
Eqs. (3.2), (3.44), (3.45) and (4.17) 

E t (k) = C lPl [e(k)/ PcS (k)] 2/3 k- 5 / 3 , (6.1) 
^)=0l±^^)- (6.2) 

It is convenient to introduce the dimensionless functions of k — k L, Ej. and E%, both normalized by 
EtiL- 1 ): 

f _ Etjk) _ E p (k) 

^-s^Fr ^-^i 31 )' (6 - 3) 
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which may be written as 



El 



Next, introduce the dimensionless ratio 

R{k) 



El 



El 



f {0) o 



2/3 



(6.4) 
(6.5) 

(6.6) 



where E® ,{ = [e^] 2 / 3 k~ 5 / 3 is the density of turbulent kinetic energy and is the energy flux in the 
particle-free case, Eq. (5.6). The ratio R(k) is larger (smaller) than unity in the case of enhancement 
(suppression) of the turbulent energy by particles. 



B. Energy flux 



Our model with local in k space parametrization of the energy flux involves the parameter of the 
closure procedure C, which has to be considered as a fit parameter which may be evaluated, for 
example, by comparison with the direct numerical simulation. Generally speaking, it is expected to 
be of the order of unity. For simplification of the qualitative analysis of the effect of particles on 
the statistics of turbulence we choose usually C = 1/4, for which we can use the analytical solution 
(5.14-5.15). The effect of C on the behavior of the function e K is qualitativly described by the 
approximation s K (C) ~ l + C [e K (C — 1)]. This is illustrated by comparison of plots £ ,« fc> r C = 0.25 
and C = 1 in Figs, la and 16. 

The plots in Fig. 1 also demonstrate the expected fact that for small the effect of the particles is 
proportional to 0. This may also clearly be seen from the balance Eq. (5.11). The balance equation 
itself also reflects a less trivial fact of saturation of the effect of the particles in the limit 3> 1; the 
beginning of this saturation is clearly seen in Fig. 1. The physical reason for the saturation is that 
the main governing parameter in the problem is the ratio of the particle energy to the energy of the 
suspension but not to the energy of the carrier fluid. As an example of the quantitative description 
of the effect of saturation consider Eq. (5.14) for e 0jK in the limit k — > oo: 



Evidently: 



£0,oo — 



1 + 



4X/TT0 n U + 1 + - ^0(1 + 0) 



1 + 



C0 



2(l + o) 



for < 1 , 



(6.7) 



(6. 



1 + ^ln 

4 
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for > 1 . 



One sees in Fig. 1 that the increase in the mass loading 0, leads to the suppression of the energy 
flux for large k ( small scales). The onset of this suppression shifts to smaller k (larger scales) 
with increasing 5, Fig. 2. To understand this, we note that r p 7 K 5 k 2 ^ 3 is an important governing 
parameter in the energy budget equation. Consequently, with increasing particle response time, the 
fluid-particle friction dissipate energy in the larger scale region. As is evident from this figure, the 
main dissipation of energy occurs in the region 0.1 < 5k 2 / 3 < 10. Therefore, for 5 > 0.1 the total 
loss in the energy decreases with further increase in 5. 



29 



y 



0.9 




C=0.25, <M 



0.8 



10 



FIG. 2: Log-Log plots of analytical solution e , K for 5 = 1, 0.1, 1CT 2 , 1CT 3 (with C = 0.25 and <j> = 1). 



C. Suppression and enhancement of turbulence in the inertial interval 



As we discussed in Sec. VIA the effect of particles on the energy distribution in suspension (with 
respect to the particle free case) may be characterized by the ratio R(k), Eq. (6.6). This effect 
is twofold: the fluid-particle friction leads to suppression of the energy flux with increasing k. 
Accordingly, e K in the numerator of Eq. (6.6) decreases towards larger n. On the other hand, for 
larger k less particles are involved in the motion and the effective density, p K , decreases with n. The 
factor = 1 in the limit 7?e s — > oo. Therefore in the inertial interval of scales the behavior of R(k) 
is defined by the strongest dependence either of e K or of p K . As we discussed, for C < 0.25 it is 
sufficient to use the analytical solutions for e 0>K , Eq. (5.14), and 7 0iK , Eq. (5.15). This gives: 



Fig. 3a demonstrates how the ratio R(k) depends on the fit parameter of our model, C, which 
appeared in the budget Eq. (5.11) in the front of the term, responsible for the fluid-particle friction. 
Clearly, the relative importance of the fluid-particle friction (with respect to the effect of the density 
variation) increases with the value of C. In particular, for C = 1 the friction dominates and R(n) < 1, 
for C = 0.25 the density variation dominates and R(k) > 1. For 0.25 < C < 1, the density of energy 
of the carrier fluid is suppressed for smaller k and enhanced towards larger k. As is clearly seen in 
the figure, the function R(k) has a minimum around some critical value n CT which depends on C 
and k. For C ~ 1 the value of k ct agrees with the expected estimate r P 7(/c cr ) = <57 Kcr ~ 1. With 
decreasing S the position of the crossover ( and of the minimum) is shifted towards larger k. It is 
evident that for k < k cv the effect of the fluid-particle friction wins, while for k > k ct the effect of 
the decrease in the effective density of the suspension is stronger. For k 3> n cr the function R(k) 
saturates. 

Notice that for C < 0.5 the analytical prediction (dashed lines) is pretty close to the "exact" 
numerical result (solid lines) indicating the qualitative validity of the analytical description of the 
effect of particles on the energy distribution in suspensions (within the model limitations). Therefore, 




(6.9) 
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FIG. 3: Log-log plots of the analytical prediction Ro(k), Eq. (6.9), (dashed lines) and the numerical result, R(k), for <fi = 1 
(remind, k = kL). Panel a: S = 0.1 and values of C indicating corresponding lines. Panel b: C = 0.5 and values of 8 
corresponding lines. 



one can find the limiting value of Roo = R{k — > oo) from Eq. (6.9): 

!> (1 + 2 6)1 2/3 




(6.10) 



5+1 + 0+ y/ffi. + 



5+1+0- y/<f>(l+<l 

The analysis of this equation shows that the largest possible enhancement of the turbulent energy 
in the inertial interval is achieved for 5 ~ 1 and increases with 0. 



D. Turbulence modification for finite Tie. 



In the previous section we discussed the mechanism of turbulent enhancement in the inertial inter- 
val caused by the density variation in the energy cascade processes. There is one more mechanism 
of the turbulent enhancement, near the viscous subrange, that may be even more important at 
moderate Tie. This effect is due to the renormalization of the kinematic viscosity in suspensions, 
v => v e fi(k), Eq. (3.4), caused again by the density variation. Since p e s(k) near the viscous cutoff, 
k ~ I/77, is larger than pf ( and consequently v e &{l/rj) < v), the extent of the inertial interval in 
suspension is therefore larger than that in the particle free case for the same energy pumping to the 
system. Within our model this effect may be described for very small particles with a response time 
smaller than the turnover time of 77-eddies. In this case the effective density is /c-independent in the 



inertial subrange, p K ■ 
Thus Eq. (6.6) yields 



e K and e 



(o) 



(for the particle free case) are given by Eqs. (5.5) and (5.6). 



Ci 



Ci 



473e f 



(6.11) 



Plots of R(k) for different T^ef are shown in Fig. 4a by dashed lines together with the numerical 
results for a quite small 5 = 1CT 3 , solid lines. The numerical results for 5 = 0.01 and S — 0.1 are 
shown in Figs. 46 and 4c. With Tfef growing above 10 6 — 10 8 we return back to the situation in the 
inertial interval, described above, see Fig. 3. For the comparison we show the plots for 73e f — > 00 in 
Fig. 4. 
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FIG. 4: The plots of the numerical results for R(k) (k = kL) for various Tiet (solid lines) for = 1 and C = Ci = 0.5. Panel 
a: 8 = 1CP 3 , dashed lines show analytical prediction, Eq. (6.11) for 8 = 0. Panels b and c: 8 = 0.01 and <5 = 0.1. 



For very small 5 the effect of particles on the turbulent statistics in the inertial interval is negligible; 
as an illustration see Fig. 4a for 5 = 1CT 3 . In this case there is only the viscous range enhancement. 
Clearly, with decreasing 7^e f this effect is more pronounced. For the moderate values of 5, there 
is a turbulence suppression in the beginning of the inertial interval, which turns into a turbulence 
enhancement in the bulk of the inertial interval, see, e.g. Fig. 36 for 5 = 0.01 and the line marked 
Tie? = oo in Fig. 46. One sees that already for Ttef = 10 4 the energy enhancement increases. This 
enlargement becomes more and more pronounced for even smaller 73ef. For 73ef = 10 2 the turbulence 
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suppression in the beginning of the inertial interval is negligible. Further development of these 
tendencies is illustrated in Fig. 4c for 5 — 0.1 



E. Brief comparison with DNS 

In order to get an analytical description of the main physical mechanisms of the particle effect on 
turbulence we used in this paper as simple as possible approximations, which nevertheless preserve 
the basic physics of the problem. In particular, we have used the differential approximation of 
the energy flux term, Eq. (4.11) with local in /c-space closure procedure, which gives a reasonable 
approximation in the extended inertial intervals of several decades. However, in the direct numerical 
simulations of turbulence in suspensions, e.g. in Ref. [4], there is almost no inertial interval, definitely 
smaller than one decade. Therefore the detailed comparison of our simple theoretical picture with 
DNS may be only qualitative. 

For such a comparison with DNS by Boivin, Simonin and Squires [4]) we re-plotted in Fig. 5 their 
Fig. 5b for the kinetic energy spectra E s (k, 0) of suspensions in the log-log coordinates (solid lines). 
The solid line, labelled by = 0, describes the particle free case, in which the energy spectrum in the 
inertial interval should be scale invariant. The K41 dependence is shown in Fig. 5 by a dash-dotted 
line, labelled by /t~ 5 / 3 . One sees that only the first half of the decade may be considered as the 
inertial interval. The viscous corrections to this dependence may be accounted for with the help of 
Eq. (5.6). Using also Eqs. (5.9) and (6.4), one gets: 



El = k-^ 



l + i(l-- 4/3 ) 



(6.12) 



With an appropriate value of TZe^ this equation reasonably approximates the numerical data almost 
in the whole decade of k, in which decays more than three orders of magnitude, see dashed line 
= 0. The chosen value T&f = 40 agrees with parameters given in Ref. [4] with an acceptable value 
of the closure parameter Ci, which enters in the definition (5.9) for T&f. With C = 13 the numerical 
solutions of Eq. (5.3) approximate well all the DNS energy spectra E { K ((p) with = 0.2, 0.5, and 1 
in the region, bounded from above by some value of k referred to as K max . In this region the spectra 
decrease from unity (at k — 1) to some values, smaller than 10~ 3 . The value of /t max decreases from 
K max ~ 14 for = spectrum to K max ~ 7 for the spectrum with 0=1. 

For k > /t max the solutions of Eq. (5.3) give too small values of the turbulent energy. As already 
discussed, this is due to the differential approximation for the energy flux, which is absolutely not 
realistic in the viscous subrange. Clearly, the larger the value of 0, the more energy is dissipated by 
the fluid-particle friction, diminishing the energy flux at the end of the inertial interval. Consequently, 
the Kolmogorov microscale rj = z/ 3 / 4 /^ 1 / 3 increases. Since n max oc 1/rj, it shifts toward smaller values. 

One observes also some deviations of the DNS data and our numerical solutions in the energy 
containing region k ~ 1. This is again related to the differential approximation for the energy flux. 
To improve the description of the particle effect on the turbulent statistics a better approximation 
for the energy transfer term is required. This calls for the more elaborated closure procedures, based 
on a proper analysis of the triad interactions. 



VII. SUMMARY 



• In this paper we propose a one-fluid dynamical model of turbulently flowing dilute suspensions 
which differs from the usual Navier-Stokes equation in two aspects: 
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FIG. 5: Log-Log plots of turbulent kinetic energy spectrum E l K ((f>) taken from [4] for (f> = 0,0.2,0.5 and 1 with 8 = 1.65 (solid 
lines), and numerical solution of Eq. (5.3) for the same values of <j> and 8 with Vef = 40 and C = 13. 

1. Instead of fluid density, pt, our model involves a k dependent effective density p e s(k) which 
varies between pi (for large k) and the mean density of suspension p s = pf(l + 0) (for small k); 

2. The model equation includes an additional damping term oc 7 P (&) which describes the fluid- 
particle viscous friction. 

• Our model may be considered as a mean-field approximation in which one uses a dynamical equation 
of motion with "effective" coefficients which depend on the statistics of the resulting stochastic 
solutions. In our case p e s(k) and 7 P (fc) are determined by the eddy turnover frequency j(k) which, 
in its turn, depends on the resulting energy distribution in the system. 

• Our model is based on the same set of assumptions (applicability of the Stokes law for the 
fluid-particle friction and space homogeneity of the particle distribution) as widely used in two- fluid 
models for suspensions. We believe that the one-fluid description of turbulent suspensions contains 
the same physics as the essentially more complicated two-fluid models. Our feeling is that a possible 
minor difference in the level of accuracy between these two models is beyond a current level of 
understanding of the problem and is definitely smaller than the "absolute" accuracy of each model 
itself. 

• In order to keep the description of the problem as simple and transparent as possible we used in 
this paper a closure procedure based on the Kolmogorov-41 dimensional reasoning with an additional 
simplification — the differential form of the energy transfer term in which the energy flux e(k) is 
evaluated locally in k-space, via the spectrum £{k) taken at the same wave-number k. This allows 
us to derive the quite simple ordinary differential equation for the energy budget in the system (5.3). 

• As a reward, our budget Eq. (5.3) allows an effective analytical analysis in various important 
limiting cases, i.e.: 

1. In the particle free case, see Eq. (5.6); 

2. For the micro-particles case (5 < -Re~ 3 / 4 ), see Eq. (5.5); 

3. For the first decades of the inertial interval (in the case 5 < 1) or in the whole inertial interval 
(if C < 1/4), see Eqs. (5.14) and (5.15); 
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4. For any reasonable values of parameters at hand, see Eqs. (5.17) and (5.18), involving one- 
dimensional integration. 

In the general case the budget equation (5.3) may be easily solved numerically. 

• We derived the analytical expression (6.9) for the dimensionless ratio R (kL), which describes 
the energy suppression and enhancement in the inertial interval of scales. 

• In Sec. VI D we described the additional "viscous" mechanism of the turbulence suppression and 
enhancement, caused by the particle effect on the extent of the inertial interval: 

1. The decrease of the effective kinematic viscosity in suspensions (due to the increase in the 
effective density for small scale motions) elongate the inertial interval. 

2. The fluid-particle friction causes a decrease of the energy flux at the viscous end of the inertial 
interval and hence shorten the inertial interval. 

The winner of this competition depends mainly on the value of 5, see, e.g. Fig. 4. 

The complicated interplay of the inertial-range and the viscous-range mechanisms of the suppres- 
sion and the enhancement of the turbulent activity in suspensions is the main topic of Sects. V and 
VI. 

• Our model successfully correlates observed features of numerical simulations. These features are 
the following: 

1. For a suspension with particles with a response time much larger than the Kolmogorov time the 
main effect of the particles is suppression of the turbulence energy of fluid eddies of all sizes (at 
the same energy input as for the particle-free case). See for instance Fig. 5, where a comparison 
with the DNS-results of Boivin, Simonin and Squires [4] is shown. 

2. For a suspension with particles with a response time comparable to or smaller than the Kol- 
mogorov time the Kolmogorov length scale of the fluid eddies will decrease and the turbulence 
energy of eddies of (nearly) all sizes increases (at the same energy input as for the particle-free 
case). This result was also reported by Druzhinin [13], who carried out DNS-simulations for 
the case of micro-particles. 

3. For a suspension with particles with a response time in between the two limiting cases mentioned 
above the energy of the larger fluid eddies is suppressed whereas the energy of the smaller eddies 
is enhanced. The cross-over between suppression and enhancement depends on the ratio of the 
particle response time and the Kolmogorov time. The strength of the effect depends on the 
mass loading. This is in agreement with the DNS-results of Sundaram and Collins [12]. 

The more detailed comparison of our approach to turbulent suspensions with the physical and 
numerical experiments requires: 

1. From DNS side more detailed analysis of joint statistics of the velocity field of the particle and 
the carried fluid; 

2. From the theoretical side an application of the more advanced non-local closure procedures, 
explicitly accounting for the triad interactions. 

• An additional advantage of our one-fluid approach is that one can use standard and well developed 
closures from analytical theory of one-phase turbulence. This fact and the relative simplicity and 
physical transparency of the one-fluid model equations may essentially help in the further progress 
towards a theory of turbulent suspensions for more realistic cases with space inhomogeneities, grav- 
itational settling, etc. 
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